
program define ngreg
	version 9.0
	if replay() {
		if ("`e(cmd)'" ~= "ngreg") error 301
		ngreg_replay `0'
	}
	else _ngreg `0'
end

program define _ngreg, eclass
	syntax varlist(min=2) [if] [in], Id(varlist min=2 max=2) [Group(varlist) SYMmetric OLS LOGit NOConstant]

	display _newline

	gettoken y x : varlist

	gettoken igroup1 igroup2 : id

	opts_exclusive "`ols' `logit'"
	if missing("`logit'") local ols "ols"
	if missing("`noconstant'") {
		tempvar cons
		generate byte `cons' = 1
		local x "`x' `cons'"
	}

	marksample touse
	markout `touse' `igroup1' `igroup2' `group' `y' `x', strok

	_rmdcoll `y' `x' if `touse', noconstant
	local x "`r(varlist)'"

	tempname coef var xepex xpx xpxinv glsm
	tempname A1 A2 xacc1 xacc3 xx xxx AA cov_net AACUM XXCUM
	tempvar ehat phat idgroup
	tempfile jjj

quietly{
	preserve
	keep if `touse'

	// prepare the data
	count if (`igroup1' == `igroup2')
	if r(N) {
		noisily display as text "note: `r(N)' diagonal observations were not used in calculations"
		drop if (`igroup1' == `igroup2')
	}
	if ~missing("`symmetric'") {
		keep `touse' `igroup1' `igroup2' `group' `y' `x'
		save `jjj'
		tempname tvar
		rename `igroup1' `tvar'
		rename `igroup2' `igroup1'
		rename `tvar' `igroup2'
		append using `jjj'
	}
	if missing("`group'") {
		tempvar group
		generate byte `group' = 1
	}

   	if missing("`ols'") {
		logit `y' `x', noconstant robust
		predict `phat'                 			// Phat is exp(XB)/(1+exp(XB))
		generate double `ehat' = `y' - `phat'   // The residual
	}
	else {
		regress `y' `x', noconstant robust
		predict double `ehat', residuals
	}

	matrix `coef' = e(b)
	matrix `var'  = e(V)

  	local nobs = e(N)
	local K    = e(df_m)
 	local invN = 1/`nobs'

	foreach i of local x {
		tempvar tmp
		generate double `tmp'=`i'*`ehat'
		local ee "`ee' `tmp'"
	}

	tempvar rank merge

		sort `group' `igroup2' `igroup1'
		generate `rank'=_n
		save `jjj', replace

		foreach i of local ee {
			tempname tmp
			rename `i' `tmp'
			local e "`e' `tmp'"
		}

	sort `group' `igroup1' `igroup2'
	replace `rank'=_n

	merge `rank' using `jjj', sort _merge(`merge') keep(`ee')

	matrix `AACUM'=J(`K',`K',0)
	matrix `XXCUM'=J(`K',`K',0)

	tabulate `group'
	local totgroup=r(r)
	egen `idgroup'=group(`group')

	forvalues ii=1/`totgroup' {
		local c =cond(mod(`ii',50)==0, "50", "_c")
		noisily display as text "." `c'

		count if `idgroup'==`ii'

		// count the number of observations in the group
		local ngroup=(1+sqrt(1+4*`r(N)'))/2

		if mod(`ngroup',1) {
			noisily display as error _n "network matrix for the group " `ii' " is not square (possibly because of missing observations)"
			exit 503
		}

		// must adjust this so that only within group
		matrix accum `xepex'=`e' if `idgroup'==`ii', noconstant
		matrix accum `xpx'  =`x' if `idgroup'==`ii', noconstant

		// first the direct correlations
		matrix define `glsm' =J(`ngroup',`ngroup',1)

		sort `igroup1' `igroup2'
		matrix glsaccum `xx'=`e' `ee' if `idgroup'==`ii', group(`igroup1') row(`igroup2') glsmat(`glsm') noconstant

		matrix `A1'=`xx'[1..`K',1..`K']
		matrix `A2'=`xx'[`K'+1...,`K'+1...]
	 noi
		matrix `xacc1'=`xx'[1..`K',`K'+1...]
  		matrix `xacc1'=0.5*(`xacc1'+`xacc1'')

		matrix accum `xxx'=`e' `ee' if `idgroup'==`ii', noconstant
		matrix `xacc3'=`xxx'[1..`K',`K'+1...]

		// now put all the pieces together
		matrix `AA'=`A1'+`A2'-`xepex'+2*`xacc1'-`xacc3'

		matrix `AACUM'=`AACUM'+`invN'*`AA'
		matrix `XXCUM'=`XXCUM'+`xpx'

	} // forvalues

	if missing("`ols'") {
		tempvar sqdphat
		generate double `sqdphat' = sqrt(`phat'*(1-`phat'))

		foreach i of local x {
        	tempvar tmp
        	generate double `tmp'=`i'*`sqdphat'
			local stf "`stf' `tmp'"
        }
		tempname dg
        matrix accum `dg'=`stf', noconstant
		matrix 		 `dg'= `invN'*`dg'

		matrix `cov_net' = `invN'*inv(`dg'*inv(`AACUM')*`dg'')
	}
	else {
		matrix `xpxinv'=`nobs'*syminv(`XXCUM')
		matrix `cov_net'=1/(`nobs'-`K')*`xpxinv'*`AACUM'*`xpxinv'
	}
} // quietly
	restore

	local x : subinstr local x "`cons'" "_cons"
	matrix rownames `cov_net' = `x'
	matrix colnames `cov_net' = `x'
	matrix rownames `var'     = `x'
	matrix colnames `var'     = `x'
	matrix colnames `coef'    = `x'

	ereturn post `coef' `var', depname("`y'") esample(`touse') obs(`nobs')

	ereturn scalar groups = `totgroup'

	local method =cond(missing("`ols'"), "logistic", "linear")

	ereturn local title "Grouped Dyadic `method' regression"
	ereturn local cmd "ngreg"
	ereturn local model "`ols'`logit'"
	ereturn local vcetype "Robust"
	ereturn local depvar "`y'"
	
	matrix `var'=e(V)
   	ereturn matrix white_var=`var'
	ereturn repost V=`cov_net'

   	ngreg_replay `0'
end

program define ngreg_replay

	display _newline
	display as text e(title) _c
  	display as text _col(40) "Number of obs"    _col(57) "= " as res %7.0f e(N)
  	display as text _col(40) "Number of groups" _col(57) "= " as res %7.0f e(groups)
  	display ""

	 tempname Tab
	.`Tab' = ._tab.new, col(5) lmargin(0) ignore(.b)
	// column        1      2     3     4     5
	.`Tab'.width	13    |12    14    16     9
	.`Tab'.strfmt    .   %11s     .     .     .
	.`Tab'.pad       .      2     2     6     1
	.`Tab'.numfmt    .  %9.0g %9.0g %9.0g %8.2f
	.`Tab'.strcolor  . result    .     .     .

    .`Tab'.sep, top
	.`Tab'.titlefmt  .      .  %10s      .     .
	.`Tab'.titles   "" "" "Robust" "" ""
	.`Tab'.titlefmt  .   %11s  %12s  %16s   %7s
	.`Tab'.titles	"`e(depvar)'" "Coef." "Std. Err." "Gr.Dyadic s.e." "t"
	.`Tab'.sep

	tempname b se se1 t cov_net white_var
*	matrix `cov_net' = e(cov_net)
	matrix `cov_net' = e(V)
	matrix `white_var' = e(white_var)

	foreach i in `: colnames e(b)' {
		scalar `b'=_b[`i']
*		scalar `se'=_se[`i']
	   	matrix `se'=`white_var'["`i'","`i'"]
		scalar `se'=`se'[1,1]
		scalar `se'=sqrt(`se')

	   	matrix `se1'=`cov_net'["`i'","`i'"]
		scalar `se1'=`se1'[1,1]
		scalar `se1'=sqrt(`se1')
		scalar `t'=`b'/`se1'
		local na = abbrev("`i'",12)

		.`Tab'.row "`na'" `b' `se' `se1' `t'
	}
	.`Tab'.sep, bottom
end
